package fem2.material.cohesive_law;

import fem2.Element;
import fem2.Variable;
import fem2.enu.VariableName;
import inf.jlinalg.Array2DMatrix;
import inf.jlinalg.ArrayVector;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;
import math2.MathUtilities;
import math2.MutableDouble;

public abstract class Mode1JirasekCohesiveLaw extends CohesiveLaw {

	private double Dnn;
	private double Dmm1;
	private double Dmm2;
	private int dim;

	public Mode1JirasekCohesiveLaw(double Gf, double ft, double Dnn, double Dmm1) {
		super(Gf, ft);
		this.Dnn = Dnn;
		this.Dmm1 = Dmm1;
		dim = 2;
	}

	public Mode1JirasekCohesiveLaw(double Gf, double ft, double Dnn, double Dmm1, double Dmm2) {
		super(Gf, ft);
		this.Dnn = Dnn;
		this.Dmm1 = Dmm1;
		this.Dmm2 = Dmm2;
		dim = 3;
	}

	/**
	 * compute value of softening law
	 * 
	 * @param x
	 *            opening norm
	 * @return
	 */
	protected abstract double computeSofteningLaw(double x);

	/**
	 * compute value of derivative of softening law
	 * 
	 * @param x
	 *            opening norm
	 * @return
	 */
	protected abstract double computeSofteningLawDerivative(double x);

	/**
	 * compute equivalent strain used in damage function
	 * 
	 * @param w
	 * @return
	 */
	private double computeEquivalentStrain(IVector w) {
		if (dim == 2) {
			double wn = w.get(0);
			double wm = w.get(1);
			return Math.sqrt(wn * wn + (Dmm1 / Dnn) * wm * wm);
		} else if (dim == 3) {
			double wn = w.get(0);
			double wm1 = w.get(1);
			double wm2 = w.get(2);
			return Math.sqrt(wn * wn + (Dmm1 / Dnn) * wm1 * wm1 + (Dmm2 / Dnn) * wm2 * wm2);
		} else {
			throw new Error("invalid dimension");
		}
	}

	/**
	 * compute gradient of equivalent strain w.r.t opening vector
	 * 
	 * @param w
	 * @return
	 */
	private double[] computeEquivalentStrainGradient(IVector w) {
		if (dim == 2) {
			double wn = w.get(0);
			double wm = w.get(1);
			double ee = computeEquivalentStrain(w);
			return new double[] { wn / ee, (Dmm1 / Dnn) * wm / ee };
		} else if (dim == 3) {
			double wn = w.get(0);
			double wm1 = w.get(1);
			double wm2 = w.get(2);
			double ee = computeEquivalentStrain(w);
			return new double[] { wn / ee, (Dmm1 / Dnn) * wm1 / ee, (Dmm2 / Dnn) * wm2 / ee };
		} else {
			throw new Error("invalid dimension");
		}
	}

	/**
	 * compute damage function
	 * 
	 * @param ee
	 *            equivalent strain
	 * @return
	 */
	private double computeDamageFunction(double ee) {
		double g = computeSofteningLaw(ee);
		return Dnn * ee / g;
	}

	/**
	 * compute derivative of damage function
	 * 
	 * @param ee
	 *            equivalent strain
	 * @return
	 */
	private double computeDamageFunctionDerivative(double ee) {
		double g = computeSofteningLaw(ee);
		double dg = computeSofteningLawDerivative(ee);
		return Dnn / g - Dnn * ee * dg / Math.pow(g, 2);
	}

	@Override
	public IVector computeJumpTraction(Element e, IVector w) {
		double ee = computeEquivalentStrain(w);
		double F = computeDamageFunction(ee);

		MutableDouble hvar = getVariable(e).getScalarVariable(VariableName.GAMMA);

		/*
		 * check for zero opening
		 */
		if (ee < MathUtilities.tol12) {
			double ft = computeSofteningLaw(0);
			if (dim == 2) {
				return new ArrayVector(ft, 0.0);
			} else if (dim == 3) {
				return new ArrayVector(ft, 0.0, 0.0);
			} else {
				throw new Error("invalid dimension");
			}
		}

		/*
		 * compute traction
		 */
		IVector t = (IVector) w.clone();

		if (dim == 2) {
			double wn = w.get(0);
			double wm = w.get(1);
			/*
			 * by checking for loading / unloading condition
			 */
			double gamma_old = hvar.getValue();
			double gamma;
			if (F < gamma_old) {
				/*
				 * unloading
				 */
				gamma = gamma_old;
			} else {
				/*
				 * loading
				 */
				gamma = F;
				/*
				 * update history variable
				 */
				hvar.setValue(gamma);
			}
			t.set(0, (Dnn / gamma) * wn);
			t.set(1, (Dmm1 / gamma) * wm);

		} //
		else if (dim == 3) {
			double wn = w.get(0);
			double wm1 = w.get(1);
			double wm2 = w.get(2);
			/*
			 * by checking for loading / unloading condition
			 */
			double gamma_old = hvar.getValue();
			double gamma;
			if (F < gamma_old) {
				/*
				 * unloading
				 */
				gamma = gamma_old;
			} else {
				/*
				 * loading
				 */
				gamma = F;

				/*
				 * update history variable
				 */
				hvar.setValue(gamma);
			}
			t.set(0, (Dnn / gamma) * wn);
			t.set(1, (Dmm1 / gamma) * wm1);
			t.set(2, (Dmm2 / gamma) * wm2);
		}

		return t;
	}

	@Override
	public IMatrix computeJumpStiffness(Element e, IVector w) {
		double ee = computeEquivalentStrain(w);
		double F = computeDamageFunction(ee);

		MutableDouble hvar = getVariable(e).getScalarVariable(VariableName.GAMMA);

		/*
		 * check for zero opening
		 */
		if (ee < MathUtilities.tol12) {
			/*
			 * set w to small value
			 */
			if (dim == 2) {
				w.set(0, MathUtilities.tol6);
				w.set(1, MathUtilities.tol6);
			} else if (dim == 3) {
				w.set(0, MathUtilities.tol6);
				w.set(1, MathUtilities.tol6);
				w.set(2, MathUtilities.tol6);
			}
			/*
			 * recompute equivalent strain and damage function
			 */
			ee = computeEquivalentStrain(w);
			F = computeDamageFunction(ee);
		}

		IMatrix C = new Array2DMatrix(dim, dim);

		if (dim == 2) {
			double wn = w.get(0);
			double wm = w.get(1);

			/*
			 * checking for loading / unloading condition
			 */
			double gamma_old = hvar.getValue();
			if (F < gamma_old) {
				/*
				 * unloading
				 */
				C.set(0, 0, Dnn / gamma_old);
				C.set(1, 1, Dmm1 / gamma_old);
			} else {
				/*
				 * loading
				 */
				double gamma = F;
				double dF = computeDamageFunctionDerivative(ee);
				double[] de = computeEquivalentStrainGradient(w);
				double tmp1 = Dnn * wn / Math.pow(gamma, 2);
				double tmp2 = Dmm1 * wm / Math.pow(gamma, 2);
				double tmp3 = dF * de[0];
				double tmp4 = dF * de[1];
				C.set(0, 0, Dnn / gamma - tmp1 * tmp3);
				C.set(0, 1, -tmp1 * tmp4);
				C.set(1, 0, -tmp2 * tmp3);
				C.set(1, 1, Dmm1 / gamma - tmp2 * tmp4);
				/*
				 * update history variable
				 */
				hvar.setValue(gamma);
			}
		} //
		else if (dim == 3) {
			double wn = w.get(0);
			double wm1 = w.get(1);
			double wm2 = w.get(2);

			/*
			 * checking for loading / unloading condition
			 */
			double gamma_old = hvar.getValue();
			if (F < gamma_old) {
				/*
				 * unloading
				 */
				C.set(0, 0, Dnn / gamma_old);
				C.set(1, 1, Dmm1 / gamma_old);
				C.set(2, 2, Dmm2 / gamma_old);
			} else {
				/*
				 * loading
				 */
				double gamma = F;
				double dF = computeDamageFunctionDerivative(ee);
				double[] de = computeEquivalentStrainGradient(w);
				double tmp1 = Dnn * wn / Math.pow(gamma, 2);
				double tmp2 = Dmm1 * wm1 / Math.pow(gamma, 2);
				double tmp3 = Dmm2 * wm2 / Math.pow(gamma, 2);
				double tmp4 = dF * de[0];
				double tmp5 = dF * de[1];
				double tmp6 = dF * de[2];
				C.set(0, 0, Dnn / gamma - tmp1 * tmp4);
				C.set(0, 1, -tmp1 * tmp5);
				C.set(0, 2, -tmp1 * tmp6);
				C.set(1, 0, -tmp2 * tmp4);
				C.set(1, 1, Dmm1 / gamma - tmp2 * tmp5);
				C.set(1, 2, -tmp2 * tmp6);
				C.set(2, 0, -tmp3 * tmp4);
				C.set(2, 1, -tmp3 * tmp5);
				C.set(2, 2, Dmm2 / gamma - tmp3 * tmp6);
				/*
				 * update history variable
				 */
				hvar.setValue(gamma);
			}
		}

		return C;
	}

	@Override
	public void Initialize(Element e) {
		Variable v = new Variable(1);
		v.registerScalarVariable(VariableName.GAMMA, true);
		setVariable(e, v);
	}

	@Override
	public void InitializeSolutionStep(Element e) {
	}

	@Override
	public void InitalizeNonlinearIteration(Element e) {
	}

	@Override
	public void FinalizeNonlinearIteration(Element e) {
	}

	@Override
	public void FinalizeSolutionStep(Element e) {
		getVariable(e).save();
	}

	@Override
	public void CalculateOutputData(Element e) {
	}

}
